Anomalous diffusion and ergodicity breaking in heterogeneous diffusion processes 
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We demonstrate the non-ergodicity of a simple Markovian stochastic processes with space- 
dependent diffusion coefficient D(x). For power-law forms D(x) ~ \x\ a , this process yield anomalous 



diffusion of the form (x (i)) 



,2/(2-a) 



Interestingly, in both the sub- and superdiffusive regimes 



we observe weak ergodicity breaking: the scaling of the time averaged mean squared displacement 
5 2 remains linear and thus differs from the corresponding ensemble average (x 2 (t)). We analyze the 
non-ergodic behavior of this process in terms of the ergodicity breaking parameters and the distri- 
bution of amplitude scatter of 8 2 . This model represents an alternative approach to non-ergodic, 
anomalous diffusion that might be particularly relevant for diffusion in heterogeneous media. 

PACS numbers: 87.10.Mn,89.75.Da,87.23.Ge,05.40.-a 



Since the early systematic studies of Perrin and Nord- 
lund [l|, single particle tracking has become a routine 
method to trace individual particles in microscopic sys- 
tems but also of animal 0] and human Q motion 
patterns. In a wide range of systems and over many 
time and length scales these measurements reveal anom- 
alous diffusion with mean squared displacement (MSD) 
(x 2 (t)) ~ tP @. Subdiffusion (0 < p < 1) was found for 
the motion of tracers in living biological cells or 
of bacteria in biofilms Supcrdiffusion (p > 1) is ob- 
served for active motion in biological cells [lfj , of mussels 
12 1 , plant lice 13 1 , or higher animals and humans [!, 0] ■ 

Of particular current interest are the ergodic proper- 
ties of a stochastic process: is the information obtained 
from time averages of a single or few trajectories rep- 
resentative for the ensemble [14, 151? Of the various 



stochastic models for anomalous diffusion, diffusion on 
random fractals has been shown to be ergodic [HI . Frac- 
tional Brownian and fractional Langevin equation mo- 
tion driven by long-ran ge c orrelated Gaussian noise are 
transiently non-ergodic [17H19l |. In contrast, subdiffusive 
Scher-Montroll continuous time random walks (CTRW) 
[2oj ] with diverging characteristic waiting time exhibit 
weak ergodicity breaking 21( , and time averages remain 



random even in the long time limit 22j . Weak ergodicity 
breaking was indeed observed for the blinking dynamics 
of quantum dots [23| , protein motion in human cell walls 
0, and lipid granule diffusion in yeast cells 

Can we come up with a simple stochastic model that 
still features such intricate non-ergodic behavior? Here 
we show that the MSD of heterogeneous diffusion pro- 
cesses (HDPs) with space dependent diffusion constant 
D(x) ~ \x\ a scales like (x 2 (t)) ^_t p with p = 2/(2 - a), 
while the time averaged MSD 8 2 scales linearly in both 
sub- (a < 0) and superdiffusive (a > 0) regimes. We 
quantify the non-ergodicity in terms of the ergodicity 
breaking parameters and the amplitude scatter distribu- 



tion of 5 2 . Wc showcase the analogies and differences 
between HDPs and other non-ergodic diffusion processes. 

Physically, a space dependent diffusivity appears a nat- 
ural description for diffusion in heterogeneous systems. 
Examples include Richardson diffusion in turbulence [24 1 
as well as mcsoscopic approaches to transport in hetero- 
genous porous media 25 and on random fractals [26| . 
Recently, maps of the local cytoplasmic diffusivities in 
bacterial and eukaryotic cells showd a highly heterogen- 
eous landscape ^fy, recalling the strongly time- var yin g 
diffusion coefficients of tracers in cells, see, e.g., Ref. [28f . 

Model. We start with the Markovian Langevin equa- 
tion for the displacement x(t) of a test particle in a 
heterogeneous medium with space dependent diffusivity 
D(x), namely, dx/dt = ^2D(x)((t). Here, £(t) is white 
Gaussian noise with 5-correlation (C(i)C(O) = ~ 
and zero mean (£(£)) = 0. In the following we in terp ret 
the Langevin equation in the Stratonovich sense [29j, to 
ensure the correct limiting transition for the noise with 
infinitely short correlation times [301. Specifically, for 
D(x) we employ the power-law form [3l| 

D(x) = D \x\ a . (1) 

The dimension of Do is [Dq] = cm 2_Q sec _1 . In the 
Stratonovich interpretation, we introduce the substitu- 
tion y = f a: [2Z)(x')] _1 / 2 (ia;', where y(t) is the Wiener 



process [29(. For the initial condition P(x,0) = S(x) a 



compressed (a < 0) or stretched (a > 0) Gaussian PDF 
MI/p-i / \ x \Vp 



P(x,t) 



cxp 



{2/p) 2 D t 



(2) 



emerges (see Supplementary Material (SM) [32j), with 
p = 2/(2 - a). The ensemble averaged MSD (x 2 (t)} = 
J x 2 P(x,t)dx becomes 



(x 2 (t)) 



r(p + i/2 ) (2 

7rV2 



2p 



(3) 
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Figure 1: Time evolution of the PDF for sub- (Left) and su- 
perdiffusion (Right). We compare simulations results (full 
lines) to the analytical result (dashed) of Eq. (J5|. The small 
discrepancy is due to the rectified form of the diffusion coef- 
ficient at the origin implemented in simulations (see text). 



Thus, for a < we find subdiffusion, and superdiffusion 
for a > 0. Brownian motion with (x 2 (t)) = 2 Dot emerges 
for a = 0, and a = 1 produces ballistic motion. The 
diffusion becomes increasingly fast when a tends to 2. 

Fig.[T]comparcs the PDF @ to simulations results (see 
below). Compared to the Gaussian shape of Brownian 
motion (a = 0), we observe a depletion in regions of 
higher diffusivity. For subdiffusion, this causes the cent- 
ral dip of the PDF, while for superdiffusion probability 
is shifted towards the origin. Note that the shape of 
the PDFs is significantly different from those of CTRW 
processes. In particular, for the subdiffusive case, the 
CTRW-PDF has a pronounced cusp at the origin [Bj]. 
Curiously, the subdiffusive shape in Fig^ T] resembles the 
propagator for retarded wave motion [33 1 . 

To characterize the HDP, we calculate the position 
autocorrelation (see SM [32j) for the case t 2 > t\, 



{x{h)x{t 2 )) 



2P+ 1 r(p + i)r(f + 1) 

^ P 2J T(§ + I) 



x[A)(t 2 -ii)] (p - 1)/2 2 Fi ( 
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(4) 



with the hypergeometric function 2 Fi(z). The Brownian 
limit (x(t)x(t + t)) = 2Dot follows for a = 0, while for 
<i = t 2 we recover Eq. ([3]). The asymptotic behavior 
of expression g]), (x(t)x(t + r)) ~ T ( P -i)/2 t ( P +i)/2 for 
r/t > 1 implies that the correlations decay with r for 
subdiffusive processes whereas they increase in the case 
of superdiffusion, similar to fractional Brownian motion 
(FBM) [13 • We also obtain the correlation of two con- 
secutive increments along x(t). For r/t 1, we find the 
simple expressions for p = 1/2 and p = 2, 



[x(t)-x(t-r)][x(t- 



-x(t)}) 



-Ty/Da/t, p=\ 

(D t) 2 , p=2 



(5) 

Indeed the occurrence of antipersistence (negative cor- 
relations) for subdiffusion and persistence for superdiffu- 
sion holds for all values of p [32|, again similar to FBM. 



We present exact results for the velocity autocorrelation 
function in SM [32 1. 

To connect to single particle tracking experiments we 
now turn to the time averaged MSD of a trajectory x(t), 



T-A 



1 



{[[x(t + A)~x(t)] 2 )dt 1 (6) 



where A is the lag time and T the length of the time series 
x(t). In Eq. (JB|) we introduced the additional average over 

N~ 



individual trajectories, (<5 2 (A)) = TV" 1 J2i=i $1{ A ) 
For A<Twe find the linear A-dependence [32j 



5 2 {A) 



r(p + i/2) 

TT 1 / 2 



2;» 



D 



A 



(7) 



Remarkably, this result holds for both sub- and superdif- 
fusive regimes, and we find the exact match (8 2 (A)} = 



(A/7 1 ) 1 



! (A)). Eq. (O is our central result. Despite 



the simple nature of the Markovian HDP, we find that it 
displays weak non-ergodicity, i.e., the scaling of time and 
ensemble averages is different, (S 2 (A)) ^ (x 2 (A)). Lin- 
ear scaling of the type of Eq. (O was found for subdiffus- 
ive CTRWs [H, but also for correlated CTRWs [H. It 
contrasts the ultrawcak non-ergodicity of superdiffusivc 
CTRWs [35|]. I n particular, the dependence on the length 
T of the time series is identical to CTRW and correlated 
CTRW subdiffusion: with increasing T the effective dif- 
fusivity D c ff ~ T p_1 decreases. In the superdiffusivc re- 
gime, this D c f[ increases with time, the particle becomes 
more mobile. For HDP processes the ageing dependence 
on T arises when the particle continues to venture into 
regions of changing diffusivities, for subdiffusive CTRWs 
this is due to increasingly longer sojourn times. 

Exploring HDPs in more detail numerically, we note 
that the Langcvin equation in the Stratonovich interpret- 
ation leads to an implicit mid-point iterative scheme for 
the particle displacement. At step i + 1, we thus have 
x i+ i - Xi = sj2D([x i+ i +Xi]/2)(y i+ i - y;), where the 
increments (j/i+i — Hi) of the Wiener process represent a 
(5-correlated Gaussian noise with unit variance. Unit time 
intervals separate consecutive iteration steps. In our nu- 
merical study we consider the two generic cases p = 1/2 
and p = 2. To avoid divergencies in the discrete scheme, 
for subdiffusion we regularize the diffusivity at x = 0, 
choosing D(x) = DoA/(A + x 2 ), such that for x 2 > A 
we recover the original relation (fT]) . For superdiffusion, to 
avoid trapping at x — and thus prohibitively expensive 
simulations, we use the shifted form D{x) — Z?o(M + 1)> 
again with the correct scaling for large 

Numerical results. From the generated trajectories, we 
compute the PDF (Fig. [TJ as well as the ensemble and 
time averaged MSDs. The MSDs are shown in Fig. [2] 
For the ensemble averaged MSD (blue curve) we observe 
excellent agreement with Eq. ([3]). In the subdiffusion 
case, the predicted behavior ([3]) is reached after less than 
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Figure 2: Ensemble and time averaged MSDs for sub- and 
superdiffusive HDPs for T — 10 5 , Do = 1, and A = 0.01. 
Expected results ([3]) and Q are shown by dashed lines. The 
time averaged MSDs are logarithmically sampled. The initial 
position is xq — 0.1 for all trajectories. For subdiffusion, the 
discrepancy of theory and simulated results for (<5 2 ) is due to 
the D(x) rectification. 



a dozen steps, such that the rectification of D(x) and 
the choice of the off-center initial position Xq are prac- 
tically negligible. For superdiffusion, an extended region 
of almost normal diffusion is observed for small particle 
displacements, as long as |x| < 1 in the rectified D{x). 
The long time behavior shows good agreement between 
theory and simulations. The time averaged MSDs and 
their mean (6 2 ) follow nicely the predicted linear scal- 
ing. Only when the lag time A approaches T the linear 
behavior levels off due to the finite trajectory length T. 

Next we study the apparent diffusion coefficients and 
scaling exponents from the individual 5 2 . The points of 
each trace S 2 were logarithmically binned. The initial 
part of the trajectory, A = 1 . . . 10 2 steps, was fitted 
by a linear law with diffusivity D\. The long-time re- 
gime, A ~ 10 2 . . . 10 4 steps, was fitted with two para- 
meters, 5 2 = DpA&. The distributions of D\, Dp, and 
/3 were then fitted with the 3-parameter gamma distri- 



bution, g(z) 



-We" 



■a I z 



, the parameters being 



fixed by the first three moments obtained from the histo- 
grams. This gamma distribution was recently proposed 
for the spread with A of S 2 of a Brownian walker (36j j . 
In addition, we extracted the amplitude scatter PDF 
<K£ = ^ 2 /(^ 2 )) °f individual S 2 around the mean (S 2 ). 

Figs. [3] and [4] show the results of this analysis, together 
with the fits to the gamma distribution, g(z). We see that 
the amplitude scatter shows a pronounced spread around 
the ergodic value £ = 1 in both sub- and superdiffusive re- 
gimes. The shape of changes only marginally with A 
in the linear region of S 2 . Broad distributions of <f>(£) are 
known from both sub- and superdiffusive CTRWs. The 
differences are, however, that for subdiffusive CTRWs 
</>(£) has a finite value at £ = due to completely im- 
mobilized particles (22|. For superdiffusive CTRWs </>(£) 
decays to zero at £ = 0, but changes significantly with A 




Figure 3: Amplitude scatter PDF for subdiffusve HDP 
(top) and distributions of the fit parameters Di, Dp, and /? 
(bottom). Parameters A and T are indicated in the panels. 
Solid envelopes are 3-parameter fits to g{z) (see text). At 
least 10 3 trajectories were analyzed for each T, with T = 10 5 . 




Figure 4: The same as in Fig. [3]but for superdiffusive HDPs. 



351 ] . The scatter distribution 0(£) is thus a good criterion 
to distinguish HDPs and CTRWs. 

More specifically, for subdiffusion the distributions of 
D\ and the amplitude scatter (/>(£) are wider, while for 
superdiffusion they show a sharper peak with a maximum 
shifted to the mean at £ = 1. The width of <f>(£) is roughly 
unity for subdiffusion and 1/2 for superdiffusion. It 
changes only slightly with T in the range T = 10 3 . . . 10 5 , 
Fig. [3] The distributions for Dp typically exhibit a long 
power- law tail with f(Dp) - Dp 12 - 1A . The distribu- 
tions for the apparent long-time exponent /3 are similar 
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Figure 5: Ergodicity breaking parameters £B (thick) and EB 
(thin blue curves) as obtained from simulations (T = 10 ). 
The theoretical asymptotes for £ B (left panel, dashed lines) 
and approximate Eq. (S9) for EB (right panel, see SM [13]). 



for sub- and superdiffusion, with (/3) p= i/2 ~ 0.86 and 
(@)p=2 ~ 0.92. Finite T-effects lead to the undershoot of 
(/3) compared to the theoretical value of unity. 

Consider now the ergodicity breaking parameter [22j 



EB = lim 

T->oo 



6 2 



5 2 



5 2 



(8) 



that provides a good measure for the dispersion of time 
averaged MSDs for different types of diffusion processes 
371 ] . The sufficient condition of ergodicity of a process 
is EB = 0. A necessary condition is that the ratio of 
the time and ensemble averaged MSD equals unity. Here 
we observe that this condition is not fulfilled for p =/= 1, 

SB = (iHAj) /(x 2 {A)) = (A/T^-p. The ergodicity 

breaking parameters EB and SB extracted from simu- 
lations data at A/T <C 1 is EB w 0.4 for sub- and 
EB ss 1.4 for superdiffusion, indicating the weakly non- 
ergodic nature of HDP processes in both regimes. The 
parameter EB follows the predicted scaling for large A 
and T values, i.e., when both ensemble and time averaged 
MSDs have converged to the theoretical scaling (Figs. [2] 
and [5} . For Brownian diffusion our simulations yield the 
correct finite-time scaling EB(A) = |— and sharp amp- 
litude scatter </>(£) ~ 5(£ — 1) at T — > oo, indicating the 
ergodic behavior for long traces and the self-averaging 
property of normal diffusion. For p = 2 at A/T — > the 
estimated theoretical value is EB = 2/3, (Eq. (S9) in SM 
(32]]). close to the value from simulations. 

The simulations show that for subfdiffusion, SB in- 
creases with A/T, while for superdiffusion it decreases 
with A/T (Fig. [5J. In both cases SB approaches unity 
as A -> T. At A/T < 1, SB for p < 1 assumes smaller 
values, indicating weaker deviations from ergodicity, in 
contrast to superdiffusion, see Fig. [5] The parameter EB 
defined via the fourth moment is overall a more robust 
characteristic of the process. For instance, it varies only 
slightly with T and A for A/T < 1 (not shown). 

Conclusions. The seemingly simple Markovian HDP 
with power-law space-dependent diffusion coefficient is 



defined in terms of a Langevin equation, which is fully 
local in both space and time. Despite this locality HDPs 
are weakly non-ergodic: the time averaged MSD <5 2 (A) 
scales linearly with the lag time A for both sub- and su- 
perdiffusion. Such a behavior was previously known for 
subdiffusion from more complex processes: CTRWs with 
diverging characteristic waitin g ti me [22] ] and CTRWs 
with correlated waiting times |34| . For superdiffusion 
only an ultraweak ergodicity breaking was reported in 
which the coefficients of the time and ensemble MSDs 
differ [35[ , while space-correlated CTRWs feature a quad- 
ratic scaling of the time averaged MSD [3~ij ]. 

The amplitudes of S 2 (A) of different trajectories show 
a broad and asymmetric distribution, another sign of 
non-ergodicity. In contrast to CTRW subdiffusion, 
however, in HDPs there is no contribution from com- 
pletely immobile trajectories. Concurrently, HDPs ex- 
hibit (anti)pcrsistcncc of the increment correlation func- 
tions, such that for subdiffusion (superdiffusion) sub- 
sequent increments preferentially have opposite (equal) 
direction. This property is similar to the noise correlator 
of FBM, moreover the velocity autocorrelator of both 
processes is strikingly similar. 

HDPs represent new tools for the description of weakly 
non-ergodic dynamics. Due to their intuitive formulation 
in terms of space-dependent diffusivities this dynamic be- 
havior of HDPs directly follows from the physical prop- 
erties of the environment. HDPs may be distinguished 
from other processes due to the difference in the amp- 
litude scatter distribution of time averaged MSDs, the 
increment autocorrelation function, as well as the two 
ergodicity breaking parameters. At the same time our 
observations for HDPs pose the interesting question for 
the minimal necessary condition for the occurrence of 
weak ergodicity breaking. 

The authors thank E. Barkai and H. Buening for dis- 
cussions and acknowledge funding from the Academy 
of Finland (FiDiPro scheme to RM) and Deutsche 
Forschungsgemeinschaft (Grant CH 707/5-1 to AGC). 



[1] 
[2] 

[3] 

[4] 

[5] 



a.cherstvy@gmail.com 



rmetzler@uni-potsdam.de 



Electronic address: 
Electronic address: 
J. Perrin, Comptes Rendus 146, 967 (1908); I. Nordlund, 
Z. Phys. Chem. 87, 40 (1914). 

C. Brauchle, D. C. Lamb, and J. Michaelis, Single 
Particle Tracking and Single Molecule Energy Transfer 
(Wiley-VCH, Weinheim, 2010); X. S. Xie et al., Annu. 
Rev. Biophys. 37, 417 (2008). 

R. Nathan et al., Proc. Natl. Acad. Sci. USA 105, 19052 
(2008); D. W. Sims et al., Nature 451, 1098 (2008). 
M. C. Gonzalez, C. A. Hidalgo, and A.-L. Barabasi, 
Nature 453, 779 (2008); D. Brockmann, L. Hufnagel, 
and T. Geisel, ibid. 439, 462 (2006). 
R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000); J. 
Phys. A 37, R161 (2004). 



5 



[6] I. Golding and E. C. Cox, Phys. Rev. Lett. 96 098102 
(2006); S. C. Weber, A. J. Spakowitz, and J. A. Theriot, 
Phys. Rev. Lett. 104, 238102 (2010); I. Bronstein et al., 
Phys. Rev. Lett. 103, 018102 (2009). 
[7] A. V. Weigel, B. Simon, M. M. Tamkun, and D. Krapf, 

Proc. Nat. Acad. Sci. USA 108, 6438 (2011). 
[8] J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber- 
Unkel, K. Berg-S0rensen, L. Oddershede, and R. Metzler, 
Phys. Rev. Lett. 106, 048103 (2011). 
[9] G. Seisenberger et al., Science 294, 1929 (2001). 
10] A. Caspi, R. Granek, and M. Elbaum, Phys. Rev. Lett. 
85, 5655 (2000); Phys. Rev. E 66, 011916 (2002). 
1] S. S. Rogers, C. v. d. Walle, and T. A. Waigh, Langmuir 

24, 13459 (2008). 
2] M. de Jager, F. J. Weissing, P. M. J. Herman, B. A. 

Nolet, and J. van de Koppel, Science 332, 1551 (2011). 
3] A. Mashanova, T. H. Oliver, V. A. A. Jansen, J. R. Soc. 

Interface 7, 199 (2010). 
4] E. Barkai, Y. Garini, and R. Metzler, Phys. Today 65(8), 
29 (2012). 

5] I. M. Sokolov, Soft Matter 8, 9043 (2012). 

6] Y. Meroz, I. Eliazar, and J. Klafter, J. Phys. A 42, 

434012 (2009); Y. Meroz, I. M. Sokolov, and J. Klafter, 

Phys. Rev. Lett. 110, 090601 (2013). 
7] W. Deng and E. Barkai, Phys. Rev. E 79, 011112 (2009). 
8] J.-H. Jeon and R. Metzler, Phys. Rev. E 85, 021147 

(2012). 

9] I. Goychuk, Phys. Rev. E 80, 046125 (2009); Adv. Chem. 

Phys. 150, 187 (2012). 
3] E. W. Montroll and G. H. Weiss, J. Math. Phys. 6, 167 

(1965); H. Scher and E. W. Montroll, Phys. Rev. B 12, 

2455 (1975). 

1] J.-P. Bouchaud, J. Phys I 2, 1705 (1992); G. Bel and E. 
Barkai, Phys. Rev. Lett. 94, 240602 (2005); A. Rebensh- 
tok and E. Barkai, ibid. 99, 210601 (2007). 

2] Y. He, S. Burov, R. Metzler and E. Barkai, Phys. Rev. 
Lett. 101, 058101 (2008); A. Lubelski, I. M. Sokolov, and 
J. Klafter, ibid. 100, 250602 (2008); S. Burov, J.-H. Jeon, 
R. Metzler, and E. Barkai, Phys. Chem. Chem. Phys. 13, 
1800 (2011). 

3] F. D. Stephani, J. P. Hoogenboom, and E. Barkai, Phys. 

Today 62(2), 34 (2009). 
4] L. F. Richardson, Proc. Roy. Soc. London Ser. A 110, 



709 (1926). 

[25] R. Haggerty and S. M. Gorelick, Water Res. Res. 31, 
2383 (1995); M. Dentz et al., Adv. Water Res. 49, 13 
(2012). 

[26] C. Loverdo et al., Phys. Rev. Lett. 102, 188101 (2009); 
B. O'Shaughnessy and I. Procaccia, Phys. Rev. Lett. 54, 
455 (1985). 

[27] B. P. English, V. Hauryliuk, A. Sanamrad, S. Tankov, N. 
H. Dekker,, and J. Elf, Proc. Natl. Acad. Sci. 108, E365 
(2011); T. Kuehn et al., PLoS One 6, e22962 (2011). 

[28] M. Platani, I. Goldberg, A. I. Lamond, and J. R. Swed- 
low, Nature Cell Biol. 4, 502 (2002). 

[29] H. Risken, The Fokker-Planck Equation, (Springer, 
Heidelberg, 1989). 

[30] See the explanation in B. J. West, A. R. Bulsara, K. 
Lindenberg, V. Seshadri, and K. E. Shuler, Physica A 
97, 211 (1979). The corresponding diffusion equation for 
the PDF P(x,t) has the symmetrized form 



dP(x,t) _ d_ 
dt dx 



[31] In our model a < 2 to ensure the growth condition for 
existence and uniqueness of the solution of a Markovian 
stochastic differential equation, see I. I. Gikhman, A. V. 
Skorokhod, Stochastic Differential Equations (Springer, 
Heidelberg, 1972) and C. W. Gardiner, Handbook of 
Stochastic Methods (Springer, Heidelberg, 2004). 

[32] Supplementary Material. 

[33] R. Metzler and J. Klafter, Europhys. Lett. 51, 492 

(2000); R. Metzler and T. F. Nonnenmacher, Phys. Rev. 

E 57 6409 (1998). 
[34] M. Magdziarz, R. Metzler, W. Szczotka, and P. 

Zebrowski, Phys. Rev. E 85, 051103 (2012); V. Tejedor 

and R. Metzler, J. Phys. A 43, 082002 (2010). 
[35] A. Codec and R. Metzler, Phys. Rev. Lett. 110, 

020603 (2013); D. Froemberg and E. Barkai, E-print 

larXiv: 1211.15391 
[36] A. Andreanov and D. S. Grebenkov, J. Stat. Mech., 

p07001 (2012). 

[37] S. M. Rytov, Yu. A. Kravtsov, and V. I. Tatarskii, Prin- 
ciples of Statistical Radiophysics 1: Elements of Random 
Process Theory (Springer, Heidelberg 1987). 



